****The dataset dataexpb is the one with meanf1 = 5 and meanf2=2 
generated for Example 2 in Wall (2009);

data a; 
infile "C:\files\dataexpb.dat";
input dummy count y1-y9 iter;
run;


proc nlmixed data = a (keep = dummy count y1-y9 iter) ; by iter;

***starting values;
*Below are the starting values from MPlus for the same model;
parms
      b10=3.5 b20=1.49 b30=.779 b40 = 1.55 b50=1.48 b60=1.17
      b11=1    b21=1    b32= 1    b42 = 1     b53=1     b63=1
      psi1 = .296 psi2 = .055 psi3=.10 psi4=.428 psi5=.884 psi6=.522 psi7=.638 psi8=.622 psi9=3.65
      gam0 = 0 gam1 = 0 gam2 = 0 gam3 = 0 mf1 = 0 mf2 = 0
      phi1=.05 phi2=.05 rho = 0 ddelta=.05;

bounds -1< rho < 1, psi1-psi9>=0, phi1-phi2>=0, ddelta>=0;

phi12=(phi1*phi2)**0.5*rho;

*****factor analysis part;

****** here is where the loading matrix is defined;
 
 f3i = gam0 + gam1*f1i + gam2*f2i + gam3*f1i*f2i+ deltai;

     mu1   = b10+ b11*f1i   ;
     mu2   = b20+ b21*f1i   ;
     mu3   = b30+ b32*f2i   ;
     mu4   = b40+ b42*f2i   ;
     mu5   = b50+ b53*f3i   ;
     mu6   = b60+ b63*f3i   ;
     mu7   =          f1i   ;
     mu8   =          f2i   ;
     mu9   =          f3i   ;

***the part of the likelihood coming from latent factor part;

ll_factpart = -.5*log(psi1) - (1/(2*psi1)) * (y1 - mu1)**2
			  -.5*log(psi2) - (1/(2*psi2)) * (y2 - mu2)**2
			  -.5*log(psi3) - (1/(2*psi3)) * (y3 - mu3)**2
			  -.5*log(psi4) - (1/(2*psi4)) * (y4 - mu4)**2
			  -.5*log(psi5) - (1/(2*psi5)) * (y5 - mu5)**2
			  -.5*log(psi6) - (1/(2*psi6)) * (y6 - mu6)**2
			  -.5*log(psi7) - (1/(2*psi7)) * (y7 - mu7)**2
			  -.5*log(psi8) - (1/(2*psi8)) * (y8 - mu8)**2
			  -.5*log(psi9) - (1/(2*psi9)) * (y9 - mu9)**2;

model dummy ~ general(ll_factpart);

random f1i f2i deltai ~ normal([mf1, mf2, 0], [phi1, 
                                               phi12, phi2, 
                                                0,     0,   ddelta]) subject = count;

run;